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Abstract 

A minimal (low-dimensional) dynamical model of the sawtooth oscillations is pre- 
sented. It is assumed that the sawtooth is triggered by a thermal instability which 
causes the plasma temperature in the central part of the plasma to drop suddenly, 
leading to the sawtooth crash. It is shown that this model possesses an isolated limit 
cycle which exhibits relaxation oscillation, in the appropriate parameter regime, 
which is the typical characteristics of sawtooth oscillations. It is further shown that 
the invariant manifold of the model is actually the slow manifold of the relaxation 
oscillation. 
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1 Introduction 



Sawtooth oscillations [1], commonly observed in current carrying, magneti- 
cally confined plasmas, are believed to be the result of resistive internal kink 
mode i.e. (m = l,n = 1) oscillation. These oscillations are characterized by 
a relatively slow rise of the electron temperature in the central region of the 
plasma column followed by a rapid drop (the crash). This is a typical nature 
of a stick-slip or relaxation oscillation [2], where the stress is slowly built up 
and then suddenly released after a certain threshold, observed in many other 
dynamical systems e.g. Portevin-Le Chatelier effect [3], a matter of interest in 
material science. 

In this work, we propose a minimal (low-dimensional) dynamical model for 
sawtooth oscillation in tokamaks, based on a transport catastrophe due to a 
thermal instability [4]. Low-dimensional or low-order dynamical system i.e. a 
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system of coupled ordinary differential equations, has several advantages over 
a detailed physical model of the actual system in the understanding of the 
global behaviour of the system. These models become powerful as they are 
supported by well developed mathematical theories which can be used to gain 
insight into the qualitative behaviour of the system such as bifurcation and 
stability [5]. Several examples of these models can be found in the context of 
understanding the behaviour of fusion plasmas ranging from the edge localized 
modes (ELM) [6] to plasma turbulence [7]. 

In spite of a great deal of experimental and theoretical research, the sawteeth 
in tokamaks continue to be a subject of exploration. For example, though nu- 
merical simulations [8] and some experimental results [9] suggest total recon- 
nection within the safety factor q — 1 surface during a sawtooth crash, there 
are also experimental evidences [10], which indicate that q remains well below 
unity during a sawtooth cycle. Apparently, there have been several important 
contributions to the understanding of the sawtooth dynamics e.g. sawteeth 
with partial reconnection, based on turbulent transport [13]. A Taylor relax- 
ation model of the sawteeth has also been considered by Gimblett and Hastie 
[14]. 

Here, we primarily focus on the dynamics of sawtooth oscillations based on 
a low-dimensional dynamical model. We rigorously prove that this dynami- 
cal model of the sawteeth based on thermal instability, besides capturing the 
important physical aspects, does exhibit well defined, isolated limit cycle oscil- 
lations, characteristics of self-excited relaxation phenomena like the sawteeth. 
Alternatively, several authors have proposed Hamiltonian models [15,16,17], 
which however can have infinite number of periodic solutions depending on 
the starting points of the evolution with the same set of physical parameters, 
which is rather inconsistent with the universal nature of sawteeth for similar 
experimental conditions. These models, despite having dissipation, possess a 
conserved quantity much like the Hamiltonian of a conjugate system. 

In Section II, we formulate the minimal dynamical model. In Section III, we 
analyse the bifurcation and stability of the system, pointing out the existence 
of an isolated and unique limit cycle and its global stability. We complete this 
section by proving the uniqueness of the limit cycle where we demonstrate 
the existence of an algebraic equation for the limit cycle. Next, in Section 
IV, we address the issue of relaxation oscillation and explore the parameter 
regime where the limit cycle exhibits sawtooth-like oscillations. In Section V, 
with the help of a renormalization group method, we prove that the invariant 
manifold of the dynamical system is indeed the slow manifold of the relaxation 
oscillation. In the Appendix, we outline the Hamiltonian approach to this 
dynamical model. 
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2 Dynamical modeling of sawtooth oscillations 



Typical sawtooth oscillations in small tokamaks (R = 1 m) exhibit linear 
growth of central temperature with few milliseconds of duration and rapid 
crash time of ~ several microseconds, in the Ohmic heating phase [18,19,10]. 
Although there are several other exotic cases viz. giant and monster sawteeth 
[11,12], we shall limit our discussion to the simpler type of sawteeth with a 
linear rise of the electron temperature. The dynamical system, controlling the 
sawtooth oscillations of the central electron temperature, then can be written 
as [4,15,20] 



3 or 

2 n ^T =E h\-^(A,T e )nT ei (1) 

~m =l{Te)A (2) 

where T e is the central electron temperature expressed in energy units, A is 
the amplitude of the oscillation, u L (A,T e ) is the rate of temperature redis- 
tribution, and 7(T e ) is the growth rate of the relevant mode. The collisional 
parallel conductivity cry oc T e 3 / 2 . In the above equations, the particle density 
n remains nearly constant during the sawtooth cycle, which is consistent with 
experimental observations. We further note that in the cases of sawtooth os- 
cillations, we are going to consider, the classical diffusion time [21] for the 
plasma current within the q < 1 volume, Tj = d e ) 2 jv ei (where r\ is the 
radius of the q = 1 surface, d e is the plasma skin depth, and v ei is the electron- 
ion collision frequency), is one order of magnitude higher than the sawtooth 
repetition time. For example, in the Alcator C-Mod (MIT) machine, typically, 
Tj ~ 80-400 msec for Ohmic regimes [22], whereas the sawtooth period (crash 
time being negligibly smaller than the period) r s ^ ~ 4 msec. Therefore it can 
be safely assumed that the current redistribution does not play any signifi- 
cant role and the corresponding E\\ remains constant. We further assume that 
the pressure redistribution parameter z/i(A,T e ) can be expressed with simple 
power laws i.e. v L {A, T e ) oc T^A a , where a and a are arbitrary constants. 

With these considerations, we note that the second term in Eq.(l) is respon- 
sible for the sawtooth crash. However, in absence of this term the general 
solution of Eq.(l) is explosive. In particular, it should include a loss term, 
which takes into account all other losses e.g. radiation loss, electron thermal 
diffusion [23], and neoclassical loss [24] etc. With a diffusive loss term, Eq.(l) 
can be rewritten as 

3 dT e 9 n m „ „ „ 3 nT e , 
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Figure 1. A typical sawtooth cycle in the Alcator C-Mod (MIT) machine (Shot No. 
#960130034 [22]). 

where te is the overall plasma energy confinement time inside the q = 1 surface 
and v\ is a proportionality constant. In fact, the presence of the diffusive loss 
term can be understood from the usual pressure gradient equation without 
taking into account the instability which causes the ejection of plasma energy 
due to sawtooth crashes, 

t? = Soh + V- X ±Vp, (4) 



where p ~ nT e is the plasma pressure, Soh is the Ohmic source term i.e. 
the heating power density per unit volume, and x± is the plasma thermal 
conductivity perpendicular to the magnetic field. Attributing all other losses 
to this diffusive term, we can write [26] 

(V • X ±V) ~ -te 1 . (5) 



Considering the fact that the plasma confinement time te scales as ~ o, 2 /xd, 
a being the minor radius of the torus and Xd, a diffusion coefficient which 
obeys Bohm-like diffusion, i.e. xd — PiXB/ a , where xb — T e /16eB, is the 
ion gyro-radius, and B is the ambient toroidal magnetic filed. Assuming that 
Tj ~ T e , the last term in Eq.(3) can be shown to scale as ~ T e 5 / 2 and Eq.(3) 
becomes 

3 -n^ = - v« L nT:A° - (3nT^, (6) 



where (3 is a constant of proportionality. Note that the above relation of elec- 
tron temperature evolution closely resembles the corresponding equation of 
sawtooth model by Kubota et. al. [27], based on transport bifurcation. 

We now consider Eq.(2) which controls the amplitude evolution of the mode. 
During one sawtooth cycle of the profile shown in Fig.l, observed in the Al- 
cator C-Mod machine at MIT [22]. Since the sawtooth crash time is much 
smaller than the sawtooth growth time, the mode amplitude can be consid- 
ered to have a 'spike-like' behavior. In the limiting case, the spike becomes a 
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singularity and can be represented by a 5 function, 



A~6(t-t ), (7) 

centered around t = t . The corresponding relation for T e can be found out 
from Eq.(l) as 

T e = t-£H(t-t ), (8) 

where £ is some constant and H(t) is the Heaviside step function. This be- 
havior of T e shows a perfect sawtooth with a crash occurring at t — to and a 
discontinuity in T e . In practice, the 'spike-like' behavior of amplitude A can 
be approximated by an exponential function, 

A~exp[-fi{t-t ) 2 ], (9) 

where \i is some large positive number. The corresponding T e is 

T e ~t-£i{v^(*-*o)}+&- (10) 

The related differential equation in A is, therefore, 

A A 

— ~-2p(t-to)A. (11) 

The coefficient of A on the right hand side of the above equation can be viewed 
as 7(T e (t)). So the complete set of dynamical equations can now be written as 



^ n ^£ = E 2 a\\ - v\nT«A a - (3nT^ 2 , (12) 

where T s is the threshold temperature for the onset of the relevant instability. 
We would like to emphasize that the model proposed by Haas and Thyagaraja 
[24], based on turbulent transport essentially reduces to the above equations 
when the neoclassical losses are discarded (they have a exclusive term to take 
care of the neoclassical loss). 

At this point, we would like to note that in reality, sawtooth oscillations in 
magnetically confined plasmas are result of evolution of a 'driven-dissipative' 
system. During a sawtooth cycle, the plasma is heated by the current through 
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Ohmic dissipation which is essentially lost to the environment after the plasma 
breakdown or the sawtooth crash. The sawtooth is repeated as the plasma 
current continues to dissipate energy into the plasma. Here, the plasma current 
is the driving mechanism of the sawtooth oscillation through which, the energy 
is being continuously dissipated. However, there is no driving source in Eq.(12) 
and the oscillations represented by Eqs. (12-13) appears to be self-sustaining 
i.e. a "self heating - dissipating - feedback - heating" cycle with a positive 
and negative damping mechanism through the terms (E\ cry - /3nT e 5 / 2 ) . The 
justification for Eq.(12) comes from the fact that we are looking only at the 
power balance equation for the confined plasma [25], not at the full burn- 
cycle of the plasma confinement. So, as far as the plasma power balance is 
concerned, we assume that there is a perpetual source of energy in the form 
of the plasma current provided externally through the loop voltage, which 
is the driving mechanism and an infinite sink which drains away the energy 
dissipated through the sawtooth crash. As the whole cycle repeats, the power 
balance equation becomes self-sustaining. 

We further note that during the sawtooth ramp phase, there is a definite 
change in the magnetic configuration, which rearranges itself after the saw- 
tooth crash either through a partial or complete reconnection [23]. As this 
magnetic reorientation causes the electron temperature to drop during the 
sawtooth crash, we have dynamically modeled this event by the second term 
{v° L nT«A a ), in Eq.(12). 



3 Stability and bifurcation — limit cycle oscillation 

In this section, we examine the possibility of existence of periodic orbits of 
Eqs. (12) and (13), in particular existence of any limit set. This is to be viewed 
in contrast to the Hamiltonian formalism (see Appendix), which always admits 
periodic solutions, though the solutions are not unique. 

We begin by expressing these dynamical equations in a generic dimensionless 
form. The electron temperature T e is normalized by the threshold temperature 
T s and define x = T e /T s . We further define the other dimensionless quantities 
as 

l^^A' = y, \-^-Tr = a. (14) 



The quantity y denotes the normalized amplitude of the mode and the quantity 
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Figure 2. Phase portrait of Eqs.(16) and (17). The arrows denote the direction of 
flow. 

Soh is the Ohmic source term corresponding to the threshold temperature T s , 



~>Oh 



3 nT s 



(15) 



the quantity cr s || being the collisional conductivity at T e = T s . With these 
definitions, the set of Eqs.(12) and (13) can be written in a generic form as 



x = x {3/2) (l-ax) -x a y, (16) 
y = p(x-l)y, (17) 

where p = <77o/£bh- Note that in writing these equations we have re-scaled 
the time as t — > So^t. As we shall see shortly that the constant a < 1 while 
p which is related to the growth rate can be quite large i.e. 1. The dots 
represent the derivatives with respect to the rescaled time t. 



In order to examine the bifurcation diagram of Eqs.(16) and (17) and existence 
of limit cycle oscillation, it is instructive to construct a phase portrait of the 
equations, which is shown in Fig. 2. The arrows on the curves show the direction 
of time. As can be seen, there are three fixed points, 'A', 'B', and 'C located 
at (x,y) = (0,0), (1, 1 — a), and (1/a, 0), respectively. Point 'A' at the origin 
is a non-isolated fixed point. Point 'C is a saddle point. The only fixed point 
of interest is the point 'B' which is either a stable or unstable point depending 
on values of a and a. We are also interested in this fixed point as this is the 
only point which is in the positive quadrant as the variables x and y may take 
only positive values. In what follows, we shall show that surrounding the fixed 
point 'B', we can construct a positively invariant region [28] lying entirely in 
the first quadrant. 
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Figure 3. A positively invariant region for Eqs.(16) and (17). 
3.1 Limit cycle oscillation 

Consider now, a set of curves fi(x) given by the equations 



fi(x) = 0, 
h{x)=C(- 



f 3 (x)=x 
U(x)=V(l-x)^ +1 , 



C(--l)+6 

a 



(18) 
(19) 

(20) 
(21) 



where C can be a large positive number, 5 is a small positive number, and T> 
and \i are positive arbitrary constants. The nullclines of the Eqs.(16) and (17) 
are given by the equations 



x — 1, 

y = x 3/2 - a (l-ax) =f(x). 



(22) 
(23) 



If we now span a region (see Fig. 3) described by the set of curves fi(x)s 
through the points (1,0), (1/a, 0), (l,C(l/a— 1)), and (x f ,y f ) — the points 
of intersection of the curves fi(x) and f(x), (x',y') being the common point 
of intersection of f3,4(x) and y = f(x), it can be shown by comparing the 
slopes of the curves fi(x)s with that of the orbit described by Eqs.(16) and 
(17), that the flow of the orbit of Eqs.(16) and (17) always crosses into the 
closed region, anticlockwise [28]. Thus the flow of the orbit of Eqs(16) and 
(17) is always bounded and the closed region in Fig. 3 constitutes a positively 
invariant region for Eqs.(16) and (17). An example set of parameters for this 
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positively invariant region are a = 0.34, a = 1/2, C — 1, // = 1, T> — 0.923, 
and 5 = 1/2, the point (x',y') = (0.322,0.287). 

What needs to shown now is the fact that for the permitted parameter regime, 
the fixed point 'B' is an unstable point. We linearize the set of equations (16) 
and (17) around the fixed point 'B' = (1,1 — a) and the eigenvalues of the 
linear part of these equations at (1, 1 — a) are 

Ai l2 = \(r± v^4A) , (24) 



where 



3 5 

T= 2 ~ 2 a ~ a ^ ^ 
A = p(l-a), (26) 

and the stability of the equilibrium point 'B' depends on the values of r and 
A. As p ^> 1 and a < 1, we can find that the point 'B' is an unstable spiral, if 

/3/2-cA 3 , , 

a < — — = a c , a < - = a c . 27 

\5/2- a) 2 y ' 



The existence of the limit cycle is now evident according to the Poincare- 
Bendixson condition [2]. 

The point 'B' also becomes unstable for a > 5/2 provided a > 1, however 
the latter is excluded by our physical conditions. As the growth rate of the 
instability should be large enough to cause a sawtooth crash, p ^> 1. This also 
ensures that the fixed point 'B' is an unstable spiral. Also, we observe that 
as a — > 1, a — > oo. In fact, at a = a c (for a certain a < a c ), a supercritical 
Hopf bifurcation occurs [2,28]. The bifurcation surface for Eqs.(16) and (17) is 
shown in Fig. 4. Note the steepening of the surface f(x, a), given by Eq.(23), as 
a — > 3/2 and the disappearance of the fixed point 'A'= (0, 0) at the bifurcation 
point. In Fig. 5, we show the limit cycle for a set of parameters a — 1, a — 
0.33, and p = 1. At the birth of the limit cycle after the supercritical Hopf 
bifurcation, the period of the limit cycle is given by the Hopf period, 

4tT 27T , . 

THopf = . | » 28 

V4A - r 2 y/ p (l - a) 
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Figure 4. Bifurcation diagram for the set of equations (16) and (17). 
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Figure 5. The limit cycle for Eqs.(16) and (17) subject to condition (27). 
3.2 Local and global stability of the limit cycle 

Though the phenomena of sawtooth oscillations in magnetically confined plas- 
mas suggests existence of a 'driven-dissipation' type limit cycle oscillation, 
experimentally observed sawtooth oscillations are, in general, perturbed by 
local irregularities and fluctuations. As limit cycles are usually characterised 
by a flow field which converges to the cycle everywhere, it is important to 
classify the orbital stability of the limit cycle in order to emphasize its response 
to noise. 

The classical method for studying orbital stability of limit cycles involves 
eigenvalues of the Poincare map [29]. Several other variants also exist in liter- 
ature, for studying local stability [30,31], among which the most common one 
is the theory of Floquet multipliers or the monodromy matrix method, where 
one follows the short-time evolution of a small perturbation 5x along the limit 
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cycle. Another widely applicable method is associated with the direction of 
slowest decay [32]. Here, we employ a method based on applying a small per- 
turbation 5x(t) to the non-stationary reference solution i.e. limit cycle of the 
flow x (t) and following the directions tangential and transverse to the flow 
along the perturbed trajectory x{t) = x (t) + 5x(t) [33]. 

Consider the dynamical system 



If the perturbations Sx(t) are small enough, the evolution can be linearised 



The linearised system (30) is transformed from the fixed basis x to the rotat- 
ing, orthogonal basis x' through a unitary rotation matrix U, 5x' = U Sx, Sx = 
U^ 1 5x', so that the linearised system (30) becomes 



d 
dt 



x(t) = F[x(t)]. 



(29) 




(30) 



— 5x' = B 5x' 



(31) 



where 



B = UAU~ l + —TJ-\ 



dt 



(32) 



is the stability matrix. 



For a two- variable system like ours [Eqs. (16-17)] 



i = f(x,y), 

y = g(x,y), 



(33) 
(34) 



the stability matrix B is given by 




(35) 



where 
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A* = [f 2 f x + g 2 g y + fg(f y + g x )]/(f 2 + g 2 ), (36) 

Xtn = l(g 2 - f)(f y + g x ) + 2fg(f x - g y )]/(f + g 2 ), (37) 
Xm = 0, (38) 

Xnn = [f 2 g y + g 2 f x - fg(f y + g x )]/(f 2 + g 2 )- (39) 

are the rates of convergence of the flow from the perturbed trajectory into the 
limit cycle. Parameter X nn represents the rate of convergence normal to the 
flow and X tt represents convergence tangential to the flow, which essentially 
averages to zero over a full cycle. The tangential and normal perturbations are 
only partially coupled as X nt = 0, which means that a tangential perturbation 
remains tangential whereas a normal perturbation couples to the tangential 
perturbations (nonzero X tn ). So a perturbation is unstable or stable depending 
on whether Ay is positive or negative. The global stability measures can be 
found out from the Lyapunov exponent, which can be obtained by phase 
averaging the A^s over an orbit with period T [33], 

Ki^^fXiA&d^ (40) 



corresponding to a particular perturbation. A negative A would imply global 
stability. In practice we use a discretised version of Eq.(40), 

m 

Aij = — Kj(t + A; At), m At = T. (41) 
m k=i 



In Fig. 6, we show the stability properties of the limit cycle oscillation of 
Eqs. (16-17) for a set of parameters a — 1, a — 0.33, and p = 100, which 
produces a relaxation type oscillation (see Fig. 8). As one can see that the 
limit cycle consists of unstable and stable parts corresponding to the normal 
perturbations. The limit cycle is globally stable as indicated by the global 
stability parameter, the Lyapunov exponent A nn = —0.00431. The plot of 
A tt shows the tangential stability which is essentially due to acceleration and 
deceleration of the flow showing a relaxation mechanism. As expected, the 
average A tt pa within the computational accuracy. The coupling matrix ele- 
ment X tn represents the effect of perpendicular perturbation on the oscillator 
phase. During X tn > 0, the perpendicular perturbation leads to advancement 
of the phase of the oscillator and X tn < indicates a phase delay. 



3. 3 Equation of the limit cycle 



In this subsection, we find out an equation for the limit cycle enclosing the 
critical point 'B' in terms of a series expansion. Here, we employ a method due 
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Figure 6. Local and global stability of the limit cycle of Eqs. (16-17) away from the 
Hopf bifurcation point. Clockwise from top, in the first plot, the unstable (dotted) 
and stable (solid) parts of the limit cycle in the phase plane are shown. The variations 
of X nn , Xtn, and \u are shown over one phase of the limit cycle in the second, third, 
and fourth plots. 

to Giacomini and Viano [34], which is based on the fact that if we consider 
the partial differential equation 



p dV_ + ^dy_ _ (dP_ + dQ\ 
dx dy \dx dx I 



(42) 



corresponding to the two-dimensional dynamical system 
x = P(x,y), y = Q(x,y), 



(43) 



there exists a unique convergent power series solution of Eq.(42) in some region 
91 containing a non-degenerate critical point P$ of node or focus (spiral) type 
and that this solution satisfies 



V(x,y)=0 



(44) 



on any limit cycle contained in 91 which encloses Pq. In our case, the equations 
(16) and (17) have an isolated limit cycle enclosing the critical point located 
at (1,1 — a). We have already seen that subject to the conditions (27), the 
equilibrium point 'B' is an unstable spiral. Therefore, there exists a unique 
series solution of Eq.(42). 
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As the cases, where an explicit series solution can be calculated analytically, 
are limited to polynomial forms for the functions P(x, y) and Q(x, y), we write 
the Eqs.(16,17) with a transformation of variable x — > x 2 , so that our system 
of dynamical equations become 



x = ^x 2 (l-ax 2 )-^x 2a - 1 y, (45) 
y = p[x 2 -\)y. (46) 

Without loss of any generality, the constant a can be set to unity, for which, 
the right hand sides of Eq. (45,46) become polynomials in x and y. Following 
Giacomini and Viano [34], we now look for a power series 

oo 

V(x,y) = Y,v n (x,y), (47) 

n=0 



where v n (x,y) is a homogeneous polynomial of degree n, 

n 

v n (x, y) = J2 c n -k,kx n ~ k y k . (48) 
In practice, we work on the truncated sum 

N 

V(x,y) = J2 v n(x,y), (49) 

n=0 

at order N and try to solve for the coefficients c n _k,k from set of simultaneous 
equations generated from partial derivatives of Eq. (42) evaluated at the critical 
point Po = (1,1 — a). The equation of the limit cycle is thus given by 

V(P ) = 0. (50) 



For a set of parameters a — 1, a — 0.33, and p — 1, satisfying conditions 
(27), we calculate the coefficients c n _k,k with the help of the computer algebra 
system of Waterloo Maple [35] and find that a close curve fairly coinciding 
the numerically calculated limit cycle exists for an order as low as iV = 8. 
The corresponding curves of the the limit cycle are shown in Fig. 7. Note that 
the value of p = 1 does not necessarily produce a sawtooth, we however have 
chosen this value so as to find a close curve for V(x, y) at the lowest possible 
order. The corresponding coefficients c n _k,k are tabulated in Table. 1. 

For the sake of completeness, we review the equivalent Hamiltonian approxi- 
mation of Eqs.(16) and (17), in the Appendix. 
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Table 1 

The coefficients c n _k,k for order N = 8. 
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Figure 7. The limit cycle for Eqs.(32) and (33) as determined by the Eq.(37) (open 
circles) and numerically calculated one (solid). 

4 Relaxation oscillation and sawtooth disruptions 

Relaxation oscillations are characterized by two very different time scales. 
As we can see from the phase portrait of Eqs.(16,17) in Fig. 2, relaxation 
oscillation, if any, must be away from the Hopf bifurcation point i.e. a = a c . 
So, we must push the system toward the boundary of the invariant region, 
defined by the equilibrium points 'A' and 'C, around the bifurcation point 
(1, 1— a). At the same time, it must be noted that far away from the bifurcation 
point i.e. for a <C a c , the orbit of the limit cycle will be swept away in the flow 
around the saddle point 'C. 
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Figure 8. Limit cycle, relaxation oscillation (sawtooth) of Eqs.(16) and (17). The 
parameters in this particular case are a = 1, a = 0.33, and p = 100. 




x 



Figure 9. The phase diagram of the limit cycle (the thin closed curve) corresponding 
to the sawtooth oscillations shown in Fig. 8. The fast and slow manifolds are shown 
in thick lines. 

A sawtooth trigger is caused by a sudden increase in a fast variable (y), as 
have already been pointed out in Sec. II, which remains almost close to zero 
for the rest of the sawtooth cycle. This ensures an almost linear rise of the 
sawtooth i.e. in the slow variable, the variable x, in our case. As shown in 
Fig. 8, Eqs.(16,17) do indeed exhibit relaxation oscillations of sawtooth-type 
for an appropriate set of parameters. The phase diagram of the corresponding 
limit cycle is shown Fig. 8. As expected the fast variable y remains nearly zero 
for most of the time except for a vary short time, when it attains very high 
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value, which causes the sawtooth crash. 



The physical mechanism of this relaxation oscillation can be understood in 
terms of similar oscillations in model sandpile [36]. The corresponding variable 
in case of a sandpile is 4>(x), the height of the pile at a position x. When sand 
is continuously added at the top of the pile, an instability occurs when the 
local slope of the pile becomes too large (V0 greater than a critical value), 
causing an avalanche, which consists of rearrangements on the surface of the 
pile. In case of sawtooth oscillation, energy from the source flows at the rate 
~ x 3 / 2 ~ a [see Eq. (16-17)], causing the temperature to rise steadily, which then 
triggers an 'avalanche' after the temperature rises above a critical value. 

4-1 The slow and fast manifolds 

On the limit cycle, if we start from a point x < 1 on the slow manifold (see 
Fig. 9), the damping term in Eq.(16) i.e. (x a y) remains very small as y ~ 0. As 
the constant a < 1, dx/dt remains positive and the slow variable x increases 
in time. At the same time dy/dt < and y further decreases as long as x < 1. 
However as x increases beyond unity, dy/dt becomes positive and y increases 
as ~ exp(t 2 ), which causes the damping term in Eq.(16) to dominate over 
the growth term, triggering a crash. As soon as x again falls below unity, y 
rapidly drops [see Eq.(17)] and the whole cycle repeats. We expect that in 
the parameter region where this relaxation oscillations occur, the sawtooth 
growth time (r st ) should be independent of the growth rate of the instability 
p ^> 1, whereas, the sawtooth crash time (r c ) should be oc p~ x . 

To examine further, we note that on the slow manifold we can change the scale 
by introducing a small quantity e which defines the fast time scale t' = t/e. 
As the growth rate p can be large and y ~ on the slow manifold, rescaling 
them as 

y = e7], p = 5/e, (51) 
we can write Eqs.(16,17) as 




ax) — ex a r], 



(52) 
(53) 



Taking the limit e — * 0, we obtain the slow manifold as 



x — x 



y = o, 



,3/ 2(1 



ax), 



(54) 
(55) 
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so that the sawtooth period is essentially given by the time that the variable 
spends on the slow manifold, 



f dx 

Tst " J xmi-axY (56) 



x 3 / 2 (l — ax) ' 



where x\^ are the points of intersection of the limit cycle with the nullcline 
y = f(x) [see Eq.(23)]. As the drop in electron temperature i.e. AT e <C 1 
during a sawtooth crash, so is the drop in the corresponding dynamical variable 
Ax and hence the amplitude of the sawtooth oscillation. Expanding Eq.(56) 
around the mean value of x — 1, we can further estimate the sawtooth period 
in terms of its amplitude, 

r st ~ (1 - a) Ax. (57) 



On the fast manifold, as y can attain large values, scaling y as rj/e and time t 
as t'e, we have from Eqs.(16,17), 



x' = ex 3/2 (l -ax) -x a 7], (58) 

rf = 6{x-l)rj, (59) 

where the (') denotes derivatives with respect to the fast time scale t' . We 
obtain the fast manifold by letting e — > 0, 



±=- X a y^ (60) 

y = p(x-l)y, (61) 
which is essentially given by 

y = p(ln x - x + 1) + y m , (62) 



for a — 1, y m being the maximum value oiy on the fast manifold which occurs 
at x — 1. Both these manifolds are shown in Fig. 9. The sawtooth crash time 
r c can be approximately estimated from Eq.(59) and (61), 



/ 



x[p(lnx — x + 1) + y, 

X e 2 



<lV (63) 



As expected, in the parameter space of relaxation oscillation (p ^> 1), the 
sawtooth period is essentially independent of the growth rate p while the 
crash time is inversely proportionate to it (see Fig. 10). Note that with the 
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Figure 10. As p becomes 3> 1, the sawtooth period r st becomes independent of p. 
The crash time r c however remains inversely proportionate to p, r c oc p~ v . The 
dotted lines are the fitted curves for r c and r = r st + r c with the index p = 0.6 and 
0.5, respectively. 

index p = 0.5, the total repetition time of the sawtooth r = r st + r c closely 
follows the Hopf period as given by Eq.(28) when the oscillations remain purely 
sinusoidal for p ~ 1. 



^.i? Simulations of sawteeth 

In Fig. 11, we show the reproduction two sets of sawtooth oscillations from the 
Alcator C-Mod tokamak at MIT [22], which has an effective major and minor 
radius of 0.68 m and 0.22 m. For the two experimental shots, (a) Shot No. 
#960130034 and (b) #960127012, the threshold temperatures T s are taken to 
be the average temperatures of the sawteeth, which are, 1.7 and 1.4 keV. The 
sawtooth time periods are respectively 4.63 and 3.6 milliseconds. The constant 
applied electric field is given by E\\ = 0.291 and 0.352 volt/m which correspond 
to the Ohmic source term Soh — 83.08 and 62.09 s^ 1 for the value of Coulomb 
logarithm In A ~ 17. The parallel plasma conductivity corresponding to the 
threshold temperature T s is given by 



The temperature drop AT e during the sawtooth crash are 0.18 and 0.14 KeV, 
for these two cases. The line-averaged plasma density in the Alcator C-Mod is 
between 1 x 10 20 and 3 x 10 20 m~ 3 and we assume a constant plasma density 
of 10 20 m~ 3 during these simulations. The Z e s is taken to be ~ 1. 



9.71 x 10 3 Z eff lnAT s 3/2 .si/: 



m 



(64) 
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Figure 11. Approximate representation of two sets of experimental sawtooth oscilla- 
tions. The numerical solutions of Eqs.(16) and (17) (thick solid lines) are superim- 
posed on the actual experimental data (the zig-zag lines). 

In these simulations, the parameter regime is dictated by the value of variable 
a in Eq.(16), which is determined from the experimental parameter to be 
~ 0.506 and 0.285 in SI units, for the shots (a) and (b). The values of the 
index a are chosen so as to get the desired relaxation oscillation with a given 
growth rate p. In the cases shown in Fig. 11, the values of a are taken to be 
0.42 and 1.1, respectively with the dimensionless growth rate p = 2000, in 
both cases. 



5 The invariant manifold 



In this section, we construct the invariant manifold of Eqs.(16) and (17), using 
the renormalization-group (RG) method [37] and show that it is indeed the 
slow manifold which governs the linear rise of the electron temperature in the 
sawtooth on the limit cycle. It has been shown elsewhere that the RG method 
can be applied to continue the unperturbative solutions of evolution equations 
valid locally around an arbitrary initial value t = to, smoothly to an arbitrary 
time t constructing the envelope of the perturbative solutions [38]. 

With the rescaling of the variables, we write Eqs. (16,17) as in Eqs. (45,46) 

x = x 2 (l — ax 2 ) — xy, (65) 
y = p(x 2 -l)y. (66) 

where we have rescaled the variables t — > It and p — > p/2. Without loss 
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of any generality we have assumed that a — 1. We further make a scale 
transformation 



x = ef , 1/ = 07, 



(67) 



where e is a small and arbitrary parameter. The evolution equations now read 



^^(l-ae 2 * 2 )-^, 



(68) 
(69) 



We now try to construct a perturbative solution u(t; to) of the above equations 
around some initial time t = to, 



V 



u(t;t Q ) = u (t;t Q ) + eui(t;t ) 



+e 2 u 2 {t;t ) + O{e 2 ). (70) 
The differential equations at successive orders to be solved are given as 



(d t -A)u (t;t )=0, 

(dt-A)u 1 (t;t ) = ($-Zorio)U2, 
(d t -A)u 2 (t;t ) = (2^ 1 -^ lVo )U 2 

+PZ0V0U1, 



(71) 
(72) 

(73) 



corresponding to the zeroth, first, and second orders. In the above equations, 




1° 'P. 



(74) 



with eigenvalues (— p, 0) corresponding to the eigenvectors U\ = (0, 1)* and 
U 2 = (1, 0)*. The zeroth order solution is given by 



u {t; t ) = d(t )U 2 + c 2 (t )e- pt C/i. 



(75) 



The first and second order solutions are respectively 



, \ ciitof TT ci(t )c2(t ) t 
Ul ^ tQ) = W^A) 2 ~ (d t -A) 6 U *> 



(76) 



21 



ci(t ) 2 (t-to)t/ 2 

t 



u 2 (t;t ) = J ds (2ci(* ) -c 2 (t )e" pi ) 



to 



X 



ci(* )(s - t ) 



+ -Ci(t )c 2 (to)(e^-e^°) 
P 



Finally, the RG equation du(t;t ) / dt \ to=t = yields 

c 1 {t)-ec 1 {t) ( Cl (t)-c 2 (t)e- pt )=0, 
c 2 (t)-e 2 p Cl (t)c 2 (t) = 0. 



(77) 



(78) 



(79) 
(80) 



Note that we have retained terms only up to the first order in the first of the 
above equations. On the limit cycle, we solve the above equations, Eqs. (79,80) 
for Ci >2 (t) and with t — > oo. Thus, from Eq.(70), the invariant manifold can be 
constructed as, 



£(*;*) = 



77(t; *) = 0, 



(81) 



with C12 as constants of integration. So, this invariant manifold in terms of 
variable x and y in Eqs. (16) and (17) is given by, 



e 2 / 2et\ 



(82) 



which essentially is the slow manifold given by Eqs. (54) and (55) in Sec. IV. 



6 Conclusion 



To conclude, we have proposed a minimal dynamical model for the sawtooth 
oscillations in current carrying plasmas e.g. tokamak plasmas. This model, 
based on the assumption that the sawtooth is triggered due to a thermal in- 
stability which ejects the plasma, shows relaxation oscillations on an isolated 
limit cycle. We have further shown that the invariant manifold of this model 
is indeed the slow manifold of the relaxation oscillation. The persistent be- 
haviour of the sawtooth oscillation across different tokamaks indicate that a 
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dynamical model based on limit cycle oscillation is consistent in contrast to 
the Hamiltonian models. 

We note that, the Hopf bifurcation, which is the key to the limit cycle behav- 
ior, around an unstable equilibrium point is based on a linear theory, which 
does not always predict the behaviour of the corresponding nonlinear model. 
However, in this case, numerical simulations seem to be in agreement with the 
predicted bifurcation analysis. 
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Appendix 

Hamiltonian formalism 

In this Appendix, we review the equivalent Hamiltonian approximation of 
Eqs.(16) and (17) [15,16]. As can be shown that, a Hamiltonian for Eqs.(16) 
and (17) exist only for a = and a = 3/2, and the equations become 



x = x W 2 \l-y), 
y = p(x-l)y. 



(83) 
(84) 



In this case, the Hamiltonian of the system can be expressed as 




(85) 



with the equivalent Hamilton's equations as 



dq &H 
dt dp ' 



dp 
~dl 



on 



dq 



(86) 
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The equivalent conjugate variables i.e 'position' and 'momentum' are, respec- 
tively, given by, 

q = lny, p = -2/y/E. (87) 



The Hamiltonian formalism ensures periodic orbits described by the family of 
curves described by the differential equation, 

dy = p(x-l)y 



which are closed. However, depending on the initial conditions of evolution, 
there can be several possible orbits rather than an unique one. Note that 
the Eqs. (83,84) are similar in nature to the famous Volterra's predator-prey 
system [2]. 
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